#!/bin/bash -e

function info() {
echo 'Program: hmmcopy with control'
echo Usage: `basename $0` "[-w window_size] [-g ref_genome_file] [-p out_prefix] <sample.bam> <control.bam>"
exit 65
}

while getopts  ":w:c:p:g:" opts
do
	case  $opts  in
	w)
		window=$OPTARG
		;;
	c)
		c_const=$OPTARG
		;;
	g)
        ref_genome=$OPTARG
        ;;
	p)
        out_prefix=$OPTARG
        ;;
	\?) info;;
	esac
done
shift $(($OPTIND - 1))

if [ -z "$2" ]; then info; fi

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------
#-----------------------------------------------
log=.log
if [ ! -e "$log" ]; then
	:> $log
fi

if [ -z $out_prefix ]; then
	out_prefix=1
fi

if [ -z "$window" ]; then
	window=200000
fi

if [ -z "$ref_genome" ]; then
	ref_genome=ref_genome
fi

samtools view -H $1 |grep '^@RG' >/dev/null

if test $? -eq 0
then
sample_name=`samtools view -H $1 |grep '^@RG' |perl -n -e '{if(/SM:(.+)\t/){print $1}}'`
fi

samtools view -H $2 |grep '^@RG' >/dev/null

if test $? -eq 0
then
control_name=`samtools view -H $2 |grep '^@RG' |perl -n -e '{if(/SM:(.+)\t/){print $1}}'`
fi



chrom=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y

base_ref_genome=`basename $ref_genome`

if [ -e "${ref_genome}.gc.${window}.wig" ]; then
	echo 2>>$log 1>&2
	echo ${ref_genome}.gc.${window}.wig exists. 2>>$log 1>&2
else
	echo 2>>$log 1>&2
	echo 2>>$log 1>&2
	echo gcCounter 2>>$log 1>&2
	echo window $window 2>>$log 1>&2
	echo $ref_genome 2>>$log 1>&2
	gcCounter -c $chrom -w $window $ref_genome > $ref_genome.gc.$window.wig
fi

if [ -e "${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig" ]; then
	echo ${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig exist. 2>>$log 1>&2
else
	echo 2>>$log 1>&2
	echo 2>>$log 1>&2
	echo mapCounter 2>>$log 1>&2
	mapCounter -c $chrom -w $window ${data_path}/encode/${base_ref_genome}.map35bp.bw \
	> ${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig
fi

bai1=$1.bai
bai2=$2.bai
base1=`basename $1 .bam`
base2=`basename $2 .bam`
control_dir=`dirname $2`
if [ ! -e "$bai1" ]; then
	
	echo samtools index $1 2>>$log 1>&2
	samtools index $1
fi
if [ ! -e "$bai2" ]; then
	
	echo samtools index $1 2>>$log 1>&2
	samtools index $1
fi


controlw=${control_dir}/${base2}.${window}.control.wig

if [ ! -e "$controlw" ]; then
	
	echo 2>>$log 1>&2
	echo 2>>$log 1>&2
	echo readCounter control 2>>$log 1>&2
	echo window $window 2>>$log 1>&2
	echo $2 2>>$log 1>&2
	readCounter -c $chrom -w $window $2 > $controlw 2>>$log
else
	echo 2>>$log 1>&2
	echo 2>>$log 1>&2
	echo $controlw exists 2>>$log 1>&2
fi


# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo readCounter sample 2>>$log 1>&2
# echo window $window 2>>$log 1>&2
# echo $1 2>>$log 1>&2
# readCounter -c $chrom -w $window $1 > $base1.$window.sample.wig 2>>$log





# Rscript ${tools_path}/script/hmmcopy.r --args \
# $base1.$window.sample.wig \
# $controlw \
# $ref_genome.gc.$window.wig \
# ${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig \
# 2>>$log 1>&2
:<<!

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo readCounter sample 2>>$log 1>&2
echo window $window 2>>$log 1>&2
echo $1 2>>$log 1>&2
readCounter -c $chrom -w $window $1 > $out_prefix.$window.sample.wig 2>>$log

!
# Rscript ${tools_path}/script/hmmcopy.r --args \
# $out_prefix.$window.sample.wig \
# $controlw \
# $ref_genome.gc.$window.wig \
# ${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig \
# 2>>$log 1>&2
cp ${tools_path}/matlab/rc_medf.m .

cat $out_prefix.$window.sample.wig |awk '{if(/^f/){print "-1"}else{print}}' > 1.rc.txt

matlab -nojvm -r "rc_medf; quit;"
mv 1.rc.txt 1.1.rc.txt
mv 1.rc_f.txt 1.1.rc_f.txt

paste $out_prefix.$window.sample.wig 1.1.rc_f.txt > 1.merge.wig.txt

cat 1.merge.wig.txt |awk -F "\t" '{if(/^f/){print $1}else{print $2}}' > $out_prefix.$window.sample.filt.wig

#-----------------------
cat $controlw |awk '{if(/^f/){print "-1"}else{print}}' > 1.rc.txt

matlab -nojvm -r "rc_medf; quit;"
mv 1.rc.txt 1.2.rc.txt
mv 1.rc_f.txt 1.2.rc_f.txt

paste $controlw 1.2.rc_f.txt > 2.merge.wig.txt

cat 2.merge.wig.txt |awk -F "\t" '{if(/^f/){print $1}else{print $2}}' > $out_prefix.$window.control.filt.wig

Rscript ${tools_path}/script/hmmcopy.r --args \
$out_prefix.$window.sample.filt.wig \
$out_prefix.$window.control.filt.wig \
$ref_genome.gc.$window.wig \
${data_path}/encode/${base_ref_genome}.map35bp.${window}.wig \
2>>$log 1>&2

